#include <lapacke.h>
#include <vector>
#include <iostream>
#include <cmath>

using namespace std;

class Function
{
public:
    virtual double operator()(double _x) = 0;
    virtual double diff(double _x){return 0;};
};

class ppFrom
{
private:
    vector<double> X, Y;
    vector<vector<double>> p;
public:
    ppFrom(vector<double>& X, vector<double>& Y)
	: X(X), Y(Y){};
    vector<vector<double>> solve_c(double m1, double m2)
	{
	    int n = X.size();
	    vector<double> mu, lambda, K;
	    //vector<vector<double>> p;
	    p.resize(n);
	    mu.resize(n);
	    lambda.resize(n);
	    K.resize(n);
	    vector<double> P = {0, 0, 0, 0};
	    K[0] = 1.0*(Y[1] - Y[0])/(X[1] - X[0]);
	    for (int i = 1; i < n-1; i++)
		{
		    mu[i] = (X[i] - X[i - 1])/(X[i + 1] - X[i - 1]);
		    lambda[i] = (X[i + 1] - X[i])/(X[i + 1] - X[i - 1]);
		    K[i] = (Y[i+1] - Y[i])/(X[i+1] - X[i]);
		}
	    double A[(n-2)*(n-2)] = {[0] = 2};
	    double B[(n-2)*1] = {[0] = -lambda[1]*m1};
	    B[n-3] = -mu[n-2]*m2;
	    for (int i = 0; i < n-3; i++)
	    {
		A[(n-2)*(i+1) + i+1] = 2;
		A[(n-2)*(i+1) + i] = mu[i+1];
		A[(n-2)*i + i+1] = lambda[i+2];
		B[i] = B[i] + 3*mu[i+1]*K[i+1] + 3*lambda[i+1]*K[i];
	    }
	    B[n-3] = B[n-3] + 3*mu[n-2]*K[n-2] + 3*lambda[n-2]*K[n-3];

	    lapack_int info,_m,_n,lda,ldb,nrhs;
	    _m = n-2;
	    _n = n-2;
	    nrhs = 1;
	    lda = n-2;
	    ldb = n-2;
	    info = LAPACKE_dgels(LAPACK_COL_MAJOR,'N',_m,_n,nrhs,A,lda,B,ldb);
	    
	    vector<double> m;
	    m.resize(n);
	    m[0] = m1;
	    m[n-1] = m2;
	    for (int i = 1; i < n-1; i++)
	    {
		m[i] = B[i-1];
	    }
	    
	    for (int i = 0; i < n-1; i++)
	    {
		P[0] = Y[i];
		P[1] = m[i];
		P[2] = (3*K[i] - 2*m[i]-m[i+1])/(X[i+1] - X[i]);
		P[3] = (m[i] + m[i+1] - 2*K[i])/((X[i+1] - X[i])*(X[i+1] - X[i]));
		p[i] = P;
	    }
	    return p;
	};
    vector<vector<double>> solve_s(double M1, double M2)
	{
	    int n = X.size();
	    vector<double> mu, lambda, K, M, F;
//	    vector<vector<double>> p;
	    p.resize(n);
	    mu.resize(n);
	    lambda.resize(n);
	    K.resize(n);
	    M.resize(n);
	    F.resize(n);
	    vector<double> P = {0, 0, 0, 0};
	    K[0] = (Y[1] - Y[0])/(X[1] - X[0]);
	    for (int i = 1; i < n-1; i++)
		{
		    mu[i] = (X[i] - X[i-1])/(X[i+1] - X[i-1]);
		    lambda[i] = (X[i+1] - X[i])/(X[i+1] - X[i-1]);
		    K[i] = (Y[i+1] - Y[i])/(X[i+1] - X[i]);
		    F[i] = (K[i] - K[i-1])/(X[i+1] - X[i-1]);
		}
	    double A[(n-2)*(n-2)] = {[0] = 2};
	    double B[(n-2)*1] = {[0] = -mu[1]*M1};
	    B[n-3] = -lambda[n-2]*M2;
	    for (int i = 0; i < n-3; i++)
	    {
		A[(n-2)*(i+1) + i+1] = 2;
		A[(n-2)*(i+1) + i] = lambda[i+1];
		A[(n-2)*i + i+1] = mu[i+2];
		B[i] = B[i] + 6*F[i+1];
	    }
	    B[n-3] = B[n-3] + 6*F[n-2];
	    
	    lapack_int info,_m,_n,lda,ldb,nrhs;
	    _m = n-2;
	    _n = n-2;
	    nrhs = 1;
	    lda = n-2;
	    ldb = n-2;
	    info = LAPACKE_dgels(LAPACK_COL_MAJOR,'N',_m,_n,nrhs,A,lda,B,ldb);
	    
	    M[0] = M1;
	    M[n-1] = M2;
	    for (int i = 1; i < n-1; i++)
	    {
		M[i] = B[i-1];
	    }
	    
	    for (int i = 0; i < n-1; i++)
	    {
		P[0] = Y[i];
		P[1] = K[i] - (M[i+1] + 2*M[i])*(X[i+1] - X[i])/6;
		P[2] = M[i]/2;
		P[3] = (M[i+1] - M[i])/(6*(X[i+1] - X[i]));
		p[i] = P;
	    }
	    return p;
	};
    vector<vector<double>> solve_n()
	{
	    solve_s(0, 0);
	    return p;
	};
};

class B3
{
private:
    double t1,t2,t3,t4,t5;
public:
    B3(double t1, double t2, double t3, double t4, double t5)
	: t1(t1), t2(t2), t3(t3), t4(t4), t5(t5)
	{};
    double operator()(double x)
	{
	    double y;
	    if (x <= t1)
		y = 0;
	    else if (x <= t2)
		y = (x-t1)*(x-t1)*(x-t1)/6;
	    else if (x <= t3)
		y = ((x-t1)*((x-t1)*(t3 - x) + (t4-x)*(x-t2))+(t5-x)*(x-t2)*(x-t2))/6;
	    else if (x <= t4)
		y = ((x-t1)*(t4-x)*(t4-x) + (t5-x)*((x-t2)*(t4-x)+(t5-x)*(x-t3)))/6;
	    else if (x <= t5)
		y = (t5-x)*(t5-x)*(t5-x)/6;
	    else
		y = 0;
	    return y;
	};
};

class B2
{
private:
    double t1,t2,t3,t4;
public:
    B2(double t1, double t2, double t3, double t4)
	: t1(t1), t2(t2), t3(t3), t4(t4)
	{};
    double operator()(double x)
	{
	    double y;
	    if (x <= t1)
		y = 0;
	    else if (x <= t2)
		y = (x-t1)*(x-t1)/2;
	    else if (x <= t3)
		y = ((x-t1)*(t3 - x) + (t4-x)*(x-t2))/2;
	    else if (x <= t4)
		y = (t4-x)*(t4-x)/2;
	    else
		y = 0;
	    return y;
	};
};

class Bsplines
{
private:
    vector<double> X,Y;
public:
    Bsplines(vector<double>& X, vector<double>& Y)
	: X(X), Y(Y)
    {};
    double solve_3(double m1, double m2, double x)
	{
	    int N = X.size() - 6;
	    vector<double> a;
	    double M[N*N];
	    double b[N];
	    M[0] = 2;
	    M[N*N-1] = -2;
	    b[0] = 3*Y[3] + m1;
	    b[N-1] = -3*Y[N+1] - m2;
	    for (int i = 1; i < N; i++)
	    {
		M[N*i + i] = M[N*i + i] + 4;
		M[N*i + i - 1] = 1;
		M[N*(i-1) + i] = 1;
		b[i] = b[i] + 6*Y[i+3];
	    }
	    
	    lapack_int info,_m,_n,lda,ldb,nrhs;
	    _m = N;
	    _n = N;
	    nrhs = 1;
	    lda = N;
	    ldb = N;
	    info = LAPACKE_dgels(LAPACK_COL_MAJOR,'N',_m,_n,nrhs,M,lda,b,ldb);

	    a.resize(N+2);
	    a[0] = b[1] - 2*m1;
	    a[N+1] = b[N-2] + 2*m2;
	    for (int i = 1; i < N+1; i++)
	    {
		a[i] = b[i-1];
	    }
	    double s=0;
	    for (int i = 0; i <= N+1; i++)
	    {
		B3 B3(X[i],X[i+1], X[i+2], X[i+3], X[i+4]);
		s=s+a[i]*B3(x);
	    }
	    return s;
	};
    double solve_2(double m1, double m2, double x)
	{
	    int N = X.size() - 4;
	    vector<double> a;
	    double M[(N-1)*(N-1)];
	    double b[N-1];
	    M[0] = 5;
	    M[(N-1)*(N-1)-1] = -1;
	    b[0] = 8*Y[2]-2*m1;
	    b[N-2] = -2*m2;
	    for (int i = 1; i < N; i++)
	    {
		M[N*i + i] = M[N*i + i] + 6;
		M[N*i + i - 1] = 1;
		M[N*(i-1) + i] = 1;
		b[i] = b[i] + 8*Y[i+2];
	    }
	    
	    lapack_int info,_m,_n,lda,ldb,nrhs;
	    _m = N-1;
	    _n = N-1;
	    nrhs = 1;
	    lda = N-1;
	    ldb = N-1;
	    info = LAPACKE_dgels(LAPACK_COL_MAJOR,'N',_m,_n,nrhs,M,lda,b,ldb);

	    a.resize(N+2);
	    a[0] = 2*m1 - b[0];
	    a[N] = 2*m2 - b[N-2];
	    for (int i = 1; i < N; i++)
	    {
		a[i] = b[i-1];
	    }
	    double s=0;
	    for (int i = 0; i <= N+1; i++)
	    {
		B2 B2(X[i],X[i+1], X[i+2], X[i+3]);
		s=s+a[i]*B2(x);
	    }
	    return s;
	};
};
